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ABSTRACT 


In  seismic  tomography,  the  reference  or  starting  models  are  commonly  one-dimensional  ( 1 D)  and  the  structural 
sensitivity  kernels  of  seismic  data  are  calculated  without  considering  the  finiteness  of  seismic  waves  in  both  time 
and  frequency  domains.  These  simplifications  result  in  a  theoretical  limit,  in  addition  to  that  from  the  data  coverage, 
on  the  structural  resolution  in  tomography  inversions.  In  small-scale  (such  as  the  crustal)  structural  imaging, 
however,  a  higher  resolution  is  desired,  and  the  1 D  reference  model  is  no  longer  adequate;  therefore,  accurate 
modeling  of  wave  propagation  in  three-dimensional  (3D)  reference  models  is  necessary.  To  address  this  problem, 
elaborate  numerical  methods,  such  as  the  finite-difference  and  spectral-element  methods,  have  been  adopted  to 
calculate  the  full-wave,  finite-frequency,  banana-doughnut,  structural  sensitivity  kernels  of  seismic  data. 

The  objective  of  this  work  is  to  develop  refined  local  and  regional  velocity  and  attenuation  models  for  selected  areas 
of  interest  ( AOls)  in  Eurasia.  We  have  carried  out  systematic  tests  to  validate  the  finite-frequency  sensitivity  kernels 
computed  by  a  staggered-grid,  finite-difference  method.  These  tests  result  in  corrections  in  the  calculation  of  the 
structural  sensitivity  kernels  as  well  as  an  important  finding  that  anomalies  in  S-wave  speed  have  a  significant 
contribution  to  P-wave  traveltime  perturbations  (Zhang  and  Shen,  2008).  Thus,  current  seismic  tomography 
practices,  in  which  P-wave  traveltimes  are  assumed  to  be  unrelated  to  S-wave  velocity  anomalies,  may  lead  to 
systematic  biases  in  inversion  results.  Furthermore,  the  different  components  of  the  same  arrival  at  the  same  receiver 
have  different  traveltime  and  amplitude  sensitivities  to  variations  in  the  velocity  structure  (Shen  et  al.,  2008).  We 
have  applied  our  finite-difference-based  full-wave,  finite-frequency  approach  to  obtain  a  refined  structural  model  in 
Eastern  Turkey. 

We  continue  the  application  of  our  finite-frequency  approach  in  the  refinement  of  the  seismic  velocity  and 
attenuation  models  and  focus  on  the  Iran  region.  Our  previous  full-wave  approach  has  been  based  on  the 
finite-difference  algorithm  in  a  Cartesian  coordinate  system.  However,  the  Iran  region  under  current  study  is 
relatively  large  so  that  the  Earth’s  curvature  effect  on  wave  propagation  cannot  be  ignored.  Thus,  we  developed  an 
alternative  finite-difference  algorithm  that  solves  the  seismic-wave  equation  in  the  polar-spherical  coordinate  system 
on  a  non-staggered  and  non-uniform  grid.  We  have  also  validated  this  algorithm  by  comparing  its  synthetic 
seismograms  with  those  calculated  by  normal-mode  summation.  Our  tomography  study  covers  the  region  39°E-66°E 
and  22°N-44°N,  with  a  depth  range  of  0-700  km.  Waveform  records  from  regional  events  and  stations  were 
collected.  For  initial  tomographic  iterations,  we  selected  31  larger  earthquakes  at  16  stations  belonging  to  the  Iran 
National  Seismic  Network.  Their  source  mechanisms  are  available  from  the  Harvard  CMT  solutions,  though  we  use 
the  revised  locations  determined  by  Engdahl  and  Bergman  (personal  communication).  We  used  the  regional  a  priori 
3D  geophysical  model  of  Pasyanos  et  al.  (2004)  as  our  starting  model  and  applied  the  spherical  finite-difference 
algorithm  to  compute  the  synthetic  seismograms.  Frequency-dependent  anomalies  in  phase  and  amplitude  are 
measured  by  the  cross-correlation  between  synthetic  and  recorded  seismograms.  Full-wave,  finite-frequency 
structural  sensitivity  kernels  of  these  measurements  are  calculated  by  the  spherical  finite-difference  scheme.  The 
final  3D  structural  model  for  the  Iran  region  will  be  obtained  through  a  non-linear  iterative  process  in  which  the 
inversion  results  from  each  iteration  will  serve  as  the  reference  model  for  the  next  iteration. 
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OBJECTIVES 

In  this  work  we  obtain  Finite-Frequency  Seismic  Tomography  (FFST)  velocity  models  for  Eurasia  and  use  them  as 
3D  reference  models  to  refine  crustal  and  shallow  upper  mantle  velocity  and  attenuation  models  for  focused  AOIs 
by  applying  the  fully  3D  tomography  approach. 

RESEARCH  ACCOMPLISHED 

Development  of  Non-Staggered,  Non-Uniform  Grid  Finite-Difference  Algorithm  in  Spherical  Geometry 

In  our  research  accomplished  in  the  previous  year,  we  focused  on  the  3D  crustal  structure  in  eastern  Turkey  and 
northern  Iran  and  carried  out  a  finite-frequency  tomography  inversion  based  on  a  3D  reference  crustal  model  based 
on  MENA  1.1  (Walter  et  al.  2000).  We  used  a  fourth-order  staggered-grid  finite-difference  algorithm  (Olsen  1994) 
to  compute  the  synthetic  seismograms  in  the  layered  crustal  model  with  variable  Moho  depth.  The  crustal  model  has 
a  minimum  shear-wave  speed  of  1 . 1  km/s.  In  our  finite-difference  simulations,  we  used  a  uniform  grid  spacing  of 
300  m,  resulting  in  reliable  waveforms  up  to  a  nominal  frequency  of  about  0.5  Hz. 

Like  most  implementations  of  the  finite-difference  method,  the  staggered-grid  finite-difference  algorithm  of  Olsen 
(1994)  we  used  in  a  previous  study  has  been  designed  for  rectangular  geometry  in  the  Cartesian  coordinate  system 
where  the  Earth’s  sphericity  has  been  ignored.  This  is  acceptable  when  the  modeling  area  is  small  so  that  the  paths 
of  regional  seismic  waves  are  relatively  short.  In  this  year’s  research,  however,  we  expanded  our  area  of  interest  to 
the  entire  area  of  Iran  and  the  surrounding  region.  As  a  result  of  the  much  larger  modeling  area,  the  paths  of  the 
seismic  wave  involved  are  often  longer  so  that  the  effect  of  the  Earth’s  curvature  on  seismic  wave  propagation  needs 
to  be  taken  into  account.  A  simple  approach  to  account  for  the  Earth’s  sphericity  is  to  make  a  so-called  Earth¬ 
flattening  transformation.  It  has  often  been  applied  in  asymptotic  treatment  of  wave  propagation  in  the  spherical 
Earth  (Chapman  and  Orcutt  1985).  However,  the  transformation  is  inexact,  and  the  errors  it  causes  in  waveforms 
have  not  been  quantitatively  studied. 

Here,  we  choose  to  develop  an  implementation  of  the  finite-difference  method  in  the  polar  spherical  coordinate 
system  so  that  the  propagation  of  seismic  waves  in  the  3D  structure  of  a  spherical  Earth  can  be  modeled  accurately. 
We  use  the  geocentric  spherical  coordinate  (GSC)  system  in  which  the  origin  is  located  at  the  center  of  the  Earth  and 
the  three  coordinate  variables  are  co-latitude  0,  azimuth  (p,  and  radius  r.  The  conversion  between  the  GSC  system 
and  the  geocentric  Cartesian  system  is 

Jx2+y2 

0  =  arctan(- - ), 

(p  =  arctan(y  /  x),  (1) 

r  =  Jx2+y2+z2. 

The  elastic  wave  equation  we  solve  can  be  written  in  the  velocity-stress  form: 


pv,t  =  V-cr  +  f,  (2) 

where  the  stress-rate  tensor  can  be  expressed  via  the  Hooke’s  law 

d  =  C:{i[(Vv)  +  (Vv)T]}.  (3) 

In  the  GSC  system,  these  equations  can  be  explicitly  written  as 


PV0J  =  —0-00,0  + - cr0(p,(p  +  O0r,r  +  “  (3 GOr  +  CJ00  COt  0  -  COt0), 

r  rsinO  r 

pv<p,t  —  — <J0(p,0  + - +  oy,r  +  —  (2<T0pCot#-30>  cot#),  (4) 

r  rsmQ  r 

fWrJ  =  -  <J0r,0  + - ; - (J(pr,(p  +  CTrr,r  +  “  (  2<Jn  ~  <J00  ~  C7(fxp  +  &6t  COt  6\ 

r  rsxnO  r 
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006  -  —  ( V0  cot  #  4-  — !— ■  V(p,(p  +  V,-  )  -I-  +  Vr  )  4-  Avr,r, 

r  sin#  r 

A  +  2^1 ,  ^1  A 

<Jqxp  = - (  Vf)  COt  #  4- - V<p%<p  +  Vr  )  4-  —  (  Vq  0  4-  iy  )  4-  /lVr  r , 

r  sin#  r 

<Trr  =  —  ( VO  COt  #  4- - Vy^  +  V,-  )  4-  —  (  VO.O  4-  Vr  )  4-  (A  4-  2/v)ty  r , 

r  sin#  r 

fj  1 

(J9<P  = - (  V(pye  + - V6,<P  -  V(p  cot  #), 

2r  sin# 

H  \  1 

=—  (— Vr.^4-  V<9,r--V6>), 

2  r  r 

//  1  1 

&(pr  ~  ~(  7  “  Vr,(p  +  V(p,r  Vp). 

2  rsin#  r 


The  free-surface  boundary  condition  can  also  be  expressed  in  the  GSC  system  as 


(Jo,  =  0,  CTyr  =  0,  <Jrr  =  0. 


(5) 


(6) 


One  of  the  drawbacks  of  many  finite-difference  implementations  in  modeling  seismic  wave  propagation  is  the 
requirement  that  the  grid  size  be  constant  in  all  three  spatial  directions  and  everywhere  in  the  modeling  volume. 
When  the  seismic  wave  speed  varies  over  a  wide  range,  the  efficiency  of  the  finite-difference  scheme  is  drastically 
reduced.  This  happens  when  the  S-wave  speed  is  very  low  near  the  surface,  such  as  in  a  sedimentary  basin,  where  a 
very  small  grid  size  is  needed.  If  a  uniform  grid  is  required,  then  the  same  small  grid  size  is  used  in  the  deeper  part 
of  the  model  where  the  wave  speed  is  much  higher  and  a  small  grid  size  may  not  be  needed.  To  address  this  problem 
in  order  to  increase  the  efficiency  of  the  finite-difference  algorithm,  we  developed  a  non-uniform  grid  approach  in 
which  the  coordinates  (0,  (p,  r)  of  the  uniform  grid  points  are  transformed  to  the  coordinates  of  points  on  a  non- 
uniform  grid,  denoted  by  (£,  q,  Q,  via  the  relations: 

#  =  #(£),  (p  =  40(77),  r  =  r(g).  (7) 

The  derivatives  of  these  transformation  relations  can  be  calculated  by  the  same  differencing  operator  D  used  in  the 
finite-difference  algorithm 


#,£  —  /}<£#,  (p,rj  —  D,q(P'  Kg  —  DgT. 

Therefore,  the  velocity-stress  form  of  the  elastic  wave  equations  (4)  and  (5)  become 

1  e  1  1  ^ 

pvo,t  =—  $.00004  + - rj>(poo<i>,n  +g,r<Jfr,g  4-  — (3cr^  4- 07*7  cot  0-gw  cot#), 

r  rsin#  r 

pv<p,t  =  —  Q(p£  4- - -thrown  +  +  —  (2cr^  cot  #  -  3oy  cot  #), 

r  rsin#  r 

1  e  1  1 

pvrj  =  —g,0(J a-,4  4- - —  +  $,r<Jrr,g  4-—  (2crrr  -  Goo  -<Tpq,+Gfr  cot#), 

r  rsin#  r 

(TOO  =—  (VflCOt#4-  — rj'tpV^n  4-Vr)4  4-  Vr)  4-  Ag,rvr,g, 

r  sin#  r 

O’w  =  — — — (Vtf  cot#  4-  )  +  —(SjWOA  4-  vr )  4-  Ag,rVrx, 

r  sin#  r 

(Trr  =  — (V0COt#4-—  TJ^V^rj  4-  Vr  )  +  —  4  Vr  )  4-  (A  4  2//)$> Vr,g , 

r  sin#  r 

u  1 

<Jo<p  =  —(ZfiVpt  +  ^—jn<(pVo,n  ~  v<p  cot#), 

2r  sin# 

l* .  1  r  1  x 

( 70t  =  —(—g,0Vr,t  4 -g,rV0,g - VO  ), 

2  r  r 

_  M ,  1  1  , 

°V  -  :  —  TJ,(pVr,rj  +  g,rV(p,g  V^), 

2  rsin#  r 


(8) 


(9) 


(10) 


with  the  free-surface  boundary  conditions  for  stress  and  velocity  components  (from  equations  6  and  10): 
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1,  , 

=  -{-q,GVr4  +  V0), 

r 

A  1  1 

Vr,g  = - ; - - - +  VQ  COt  6  +  ——Tl^V^rj  +  2vr  ), 

A  +  2/ur  sin# 

1  1 

= - —  n,<pVr,rj  +-V<p. 

rsin#  r 


(ID 


We  have  implemented  this  spherical,  non-uniform  grid  approach.  A  number  of  experiments  have  been  conducted  to 
validate  the  code  by  comparing  the  simulation  results  with  the  synthetics  calculated  by  normal-mode  summation. 
Figure  1  shows  a  few  examples  of  the  comparisons  of  the  finite-difference  and  normal-mode  synthetic  seismograms. 
This  alternative  implementation  of  the  finite-difference  method  enables  us  to  efficiently  and  accurately  model  the 
propagation  of  seismic  waves  in  a  spherical  Earth  with  3D  variations  in  seismic  wave  speeds.  As  a  first  application 
of  this  finite-difference  algorithm,  we  use  it  to  conduct  a  full-wave,  finite-frequency  tomography  for  the  Iran  region. 

Full-Wave,  Finite-frequency  Tomography  for  the  Iran  Region 

Iran  is  located  in  a  seismically  active  region  with  frequent  and  widely  distributed  moderate-sized  earthquakes  (see 
Figure  2).  Yet  there  has  been  a  paucity  of  detailed  investigation  of  the  seismic  structure  of  the  crust  and  uppermost 
mantle  in  this  region.  Occasional  seismic  interests  in  this  region  included  crustal  studies  using  the  Pn  waves 
(e.g.,  Chen  et  al.  1980;  Asudeh  1982;  Hearn  and  Ni  1994;  Ritzwoller  et  al.  2002)  and  uppermost  mantle  tomography 
studies  (e.g.,  Ritzwoller  and  Levshin  1998;  Ritzwoller  et  al.  1998;  Maggi  and  Priestley  2005).  In  all  the  tomography 
studies,  the  Iran  region  is  merely  a  portion  of  the  larger  model,  and  is  characterized  by  a  single  swath  of 
low-velocity  region  without  internal  structure.  Therefore,  a  more  focused  study  using  local  earthquakes  and  stations 
is  necessary  to  obtain  a  refined  3D  structure. 

We  selected  the  region  between  22°N  and  44°N  and  between  39°E  and  66°E  for  our  full-wave,  finite-frequency 
tomography.  For  the  reference  or  starting  model,  we  chose  the  regional  a  priori  3D  geophysical  model  of  Pasyanos 
et  al.  (2004).  This  3D  model  is  constructed  by  combining  the  results  from  a  number  of  published  structural  studies 
rather  than  from  tomographic  inversions.  The  resulting  model  includes  density  and  the  velocities  and  quality  factors 
of  P  and  S  waves  specified  on  a  l°xl°  horizontal  grid  and  at  24  depths  from  0  to  660  km.  Figure  3  shows  a  few  map 
views  and  profiles  of  the  3D  reference  model. 

We  used  the  non-staggered  grid  finite-difference  algorithm  described  in  the  previous  section  to  compute  the 
synthetic  seismograms  as  well  as  the  Frechet  kernels  for  tomography  inversions.  We  stripped  the  water  from  the 
reference  model  by  uniformly  extending  the  structure  right  beneath  the  water  up  to  the  surface  and  imposed  a 
minimum  shear-wave  speed  of  1.2  km/s.  Since  most  of  the  event-station  paths  are  off  the  water/sediment  covered 
areas  (Figure  3),  this  approximation  reduces  computation  with  a  minimum  loss  of  accuracy  in  synthetic  waveforms. 
A  more  refined  grid  can  be  applied  to  sediment-covered  areas  in  following  iterations  of  solution.  In  our  finite- 
difference  simulations,  we  used  a  time  step  of  0.12  s,  a  horizontal  grid  spacing  of  0.02°,  and  a  vertically  variable 
grid  spacing  of  1-7.5  km,  as  shown  in  Figure  4,  resulting  in  reliable  surface  waveforms  up  to  a  nominal  period  of 
about  15  s  for  event-station  paths  through  the  water/sediment  covered  areas  and  ~8  s  for  those  elsewhere.  For 
P  waves,  the  period  of  reliable  waveforms  extends  to  5  s. 

We  measured  the  frequency-dependent  phase  anomalies  between  the  synthetic  and  recorded  seismograms  in  several 
narrow  frequency  bands  (0.0125-0.025,  0.025-0.05,  0.05-0.1, 0. 1-0.2  Hz)  and  a  broad  frequency  band 
(0.0125-0.2  Hz)  by  waveform  cross-correlation.  For  the  same  arrival  (say,  Pn),  we  made  measurements  on  all  useful 
components  of  the  three-component  records,  since  the  different  components  of  the  same  arrival  provide  independent 
constraints  on  the  Earth  structure  (Shen  et  al.,  2008).  Figure  5  shows  an  example  of  the  band-pass-filtered  observed 
and  synthetic  Pn  and  Rayleigh  waveforms.  We  obtained  over  2,400  high-quality  measurements  for  the  various 
arrivals,  frequency  bands,  and  components. 

In  order  to  compute  the  Frechet  kernels  of  the  measurements  with  respect  to  the  model  parameters,  we  have  adopted 
the  strain  Green  tensor  (SGT)  database  approach  developed  by  Zhao  et  al.  (2005;  2006).  We  have  calculated  the 
receiver  strain  Green  tensors  (RSGTs)  for  the  16  stations  we  have  selected  in  the  INSN  network  as  well  as  the 
source  strain  Green  tensors  (SSGTs)  for  31  events  (see  Zhao  et  al.  2006  for  details  of  the  Frechet  kernels  using  the 
RSGTs  and  SSGTs).  The  numerical  calculations  were  carried  out  with  45  processors  on  the  Linux  cluster  at 
URI/GSO.  The  resulting  SGT  database  occupies  a  disk  space  of  about  13  Terabytes.  The  Frechet  kernels  for  all  the 
measurements  can  be  efficiently  calculated  using  the  SGTs  in  the  database  in  an  iterative  full-wave,  finite-frequency 
3D  tomography  for  the  Iran  region. 
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Figure  6  shows  an  example  of  the  sensitivity  kernels  for  the  Pn  and  Rayleigh  waves  in  Figure  5.  Note  the 
sensitivities  of  the  Pn  traveltime  to  S  wave  speed.  To  a  lesser  extent,  the  Rayleigh  wave  phase  is  also  affected  by  the 
P  wave  speed  near  the  surface.  Using  numerical  and  analytical  solutions,  Zhang  and  Shen  (2008)  show  that  for 
Finite-frequency  waves,  S  wave-speed  perturbations  may  have  significant  effects  on  P  waveforms.  At  short 
source-receiver  distances  and  near  the  source  and  receiver,  the  cross-dependence  of  P  waveforms  on  S  wave-speed 
is  substantial.  This  finite-frequency  effect  has  important  implications  for  seismic  tomography,  particularly  for 
regional  and  local  earthquake  tomography,  in  which  the  volume  of  interest  is  often  within  several  10  s  of  wavelength 
from  receivers  or  both  sources  and  receivers.  Even  for  regions  far  from  sources  and  receivers,  neglecting  this 
cross-dependence  means  that  the  estimates  of  the  P  wave-speed  near  the  source  and  receiver  are  likely  biased,  and 
such  errors  may  propagate  into  the  rest  of  the  model  in  tomographic  inversions. 

CONCLUSIONS  AND  RECOMMENDATIONS 

We  have  made  great  efforts  in  developing  and  validating  finite-frequency  sensitivity  kernels  for  3D  velocity  and 
attenuation  reference  models.  We  have  implemented  a  non-staggered-grid  finite-difference  algorithm  in  spherical 
geometry  which  allows  us  to  conduct  fiill-wave,  finite-frequency  non-linear  3D  tomography  in  a  large  region  where 
the  Earth’s  sphericity  effect  on  seismic  wave  propagation  cannot  be  neglected.  We  have  calculated  the  strain  Green 
tensors  for  the  INSN  stations  and  selected  earthquakes.  Using  waveform  cross-correlation,  we  measured  over  2,400 
body  and  surface  wave  phase  anomalies  relative  to  synthetic  waveforms  of  an  initial  3D  reference  model.  We  are 
applying  the  full-wave-field  kernels  to  obtain  a  refined  3D  model  for  the  Iran  region.  Finally,  substantial  efforts  are 
needed  to  further  develop  the  method  for  large-scale  3D  models  and  models  with  anisotropy. 
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Figure  1.  Comparisons  between  two  sets  of  3-component  synthetics  calculated  by  the  spherical  finite-difference  algorithm 
developed  in  this  study  and  by  normal-mode  summation.  These  are  for  a  surface  source  located  on  the  equator 
with  strike,  dip  and  rake  angles  of  45°,  45°  and  45°,  respectively,  and  surface  receivers  to  the  east  of  the  source 
at  5°  (a)  and  10°  (/>)  distances. 
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Figure  2.  Relief  map  of  the  Iran  region  under  study.  Red  stars  show  the  epicenters  of  the  31  selected  earthquakes,  and 

blue  triangles  show  the  INSN  stations  used  in  our  tomography  inversions.  Small  red  dots  mark  the  locations  of 
the  earthquakes,  for  which  we  have  collected  waveforms.  Our  modeling  volume  is  in  the  area  outlined  by  the 
rectangle  and  has  a  depth  range  of  0  to  700  km. 
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Figure  3.  Plots  of  the  3D  reference  model  (Pasvanos  et  al.  2004).  (1st  row)  Map  views  of  Vp  at  0  km  and  30  km  depths.  (2nd 
row  )  Map  views  of  Vs  at  0  km  and  30  km  depths.  (3rd  row)  Profiles  of  Vp  along  34°N  and  57°E.  (4th  row) 

Profiles  of  Vs  along  34°N  and  57°E. 


323 


2008  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


Figure  4.  Depth  variation  of  the  radial  grid  size  used  in  the  spherical  non-uniform  grid  finite-difference  simulations  in  the 
3D  reference  model. 
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Figure  5  Pn  (left  panels)  and  Rayleigh  (right  panels)  waves  from  14  May,  2005  earthquake  ((Lat.  120.773, 
Long.  56.899)  recorded  on  the  vertical  component  at  station  ASAO  (Figure  2).  The  synthetic 
waveforms  calculated  using  the  finite-difference  code  in  a  spherical  coordinate  system  (top  traces) 
are  plotted  for  comparison.  The  waveforms  have  been  band-pass  filtered  between  0.05-0.1  Hz.  The 
vertical  lines  mark  the  time  windows,  within  which  cross-correlation  measurements  and 
corresponding  sensitivity  kernels  (Figure  6)  are  calculated. 
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Figure  6  (Top  panels)  Traveltime  Sensitivity  kernels  of  the  Pn  wave  in  Figure  5  to  P  wave-speed  (left)  and  S 
wave-speed  (right)  perturbations.  The  kernels  are  so  defined  that  low-velocity  anomaly  in  the 
negative  sensitivity  region  (red  colors)  results  in  a  phase  delay,  while  the  same  anomaly  in  the 
positive  sensitivity  region  corresponds  to  an  earlier  arrival.  (Lower  panels)  Phase  sensitivity  kernels 
of  the  Rayleigh  wave  in  Fig.  5  to  P  wave-speed  (left)  and  S  wave-speed  (right)  perturbations.  Note 
the  vertical  and  color  scales  are  different  from  those  for  the  Pn  wave. 
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